OHIBC: Plot Scores

knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'Figs/',
                      echo = TRUE, message = FALSE, warning = FALSE)

library(ohicore) ### devtools::install_github('ohi-science/ohicore')

source('~/github/ohibc/src/R/common.R')

dir_ohibc  <- '~/github/ohibc'
dir_calc   <- file.path(dir_ohibc, 'calc_ohibc')
dir_master <- file.path(dir_calc, 'master')

source(file.path(dir_calc, 'calc_scores_fxns.R'))

### provenance tracking
# library(provRmd); prov_setup()

0.1 Flower plots, overall, by year

source('~/github/ohibc/src/R/plot_flower.R')
index_scores <- read_csv(file.path(dir_calc, 'scores_all.csv')) %>%
  filter(region_id == 0) %>%
  filter(dimension %in% c('score')) %>%
  filter(goal != 'Index')

plot_wts <- read_csv(file.path(dir_calc, 'conf/goals.csv')) %>%
  select(order = order_hierarchy,
         goal, parent, goal_label = name_flower, 
         weight)

index_scores <- index_scores %>%
  left_join(plot_wts, by = 'goal') %>%
  arrange(order) %>%
  filter(!goal %in% plot_wts$parent)

for(yr in index_scores$year %>% unique() %>% sort()) {
  # yr <- 2016
  
  scores_tmp <- index_scores %>%
    filter(year == yr)
  
  index_flower <- plot_flower(scores_tmp, show_plot = FALSE) +
    labs(title = paste0('OHIBC ', yr))

  print(index_flower)
  
}

0.2 Figure: Past “likely future status” compared to current “status”

scores_all <- read_csv(file.path(dir_calc, 'scores_all.csv')) %>%
  filter(region_id == 0) %>%
  spread(dimension, score) %>%
  group_by(goal, region_id) %>%
  mutate(pred_future_status = lag(future, 5)) %>%
  arrange(year)


score_compare_plot <- ggplot(scores_all, aes(x = pred_future_status, y = status, color = goal)) +
  ggtheme_plot() +
  geom_abline(slope = 1, intercept = 0, color = 'darkred') +
  geom_point(size = 2, aes(label = year)) +
  geom_path(aes(group = goal, label = year))

plotly::ggplotly(score_compare_plot)

scores_rgn <- read_csv(file.path(dir_calc, 'scores_all.csv')) %>%
  filter(region_id != 0) %>%
  spread(dimension, score) %>%
  group_by(goal, region_id) %>%
  mutate(pred_future_status = lag(future, 5)) %>%
  arrange(year) %>%
  left_join(get_rgn_names(), by = c('region_id' = 'rgn_id'))

goals <- scores_rgn$goal %>% unique()

for(goalname in goals) { # goalname <- goals[1]
  tmp_df <- scores_rgn %>%
    filter(goal == goalname)
  
  tmp_labels <- tmp_df %>%
    group_by(rgn_name, rgn_code) %>%
    filter(!is.na(pred_future_status) & !is.na(status)) %>%
    summarize(x = last(pred_future_status),
              y = last(status))
    # summarize(x = mean(pred_future_status, na.rm = TRUE),
    #           y = mean(status, na.rm = TRUE))
    

  rgn_score_compare_plot <- ggplot(tmp_df, 
                               aes(x = pred_future_status, y = status, color = rgn_name)) +
    ggtheme_plot() +
    geom_abline(slope = 1, intercept = 0, color = 'darkred') +
    geom_point(size = 2, aes(label = year)) +
    geom_path(aes(group = rgn_name, label = year), alpha = .4) +
    geom_text(data = tmp_labels, aes(x, y, label = rgn_code), color = 'grey30') +
    labs(title = goalname,
         color = goalname)
  
  print(rgn_score_compare_plot)
}

0.3 Regression? goal-by-goal

For each goal, regress the future status (+5 years) against current status, trend, pressures, and resilience; this would give an idea of coefficients in the OHI LFS model; these likely to differ for each goal. These do not apply to supragoals (since no pressures/resilience applied at this level), only standalone and subgoals.

Assume coefficients will be similar region-by-region so compare all regions simultaneously (exclude overall region)

Goals more susceptible to pressures or amenable to resilience, worth looking at?

### obs_future_status is status five years hence; use lead() to pull future
### status to current comparison year
scores_rgn <- read_csv(file.path(dir_calc, 'scores_all.csv')) %>%
  filter(region_id != 0) %>%
  spread(dimension, score) %>%
  select(goal, region_id, year, 
         pressures, resilience, status, trend) %>%
  group_by(goal, region_id) %>%
  mutate(obs_future_status = lead(status, 5)) %>%
  arrange(year) %>%
  left_join(get_rgn_names(), by = c('region_id' = 'rgn_id')) 

goals <-  c('AO',  'FIS', # 'MAR', # 'SAL', 
            'CPP', 'CSS', 'CW',
            'HAB', 'SPP', 'ICO', 'LSP', 
            # 'LEF', 'LEO', 
            'TR')

lfs_lm_all <- data.frame()
stat_tr_all <- data.frame()

for(goalname in goals) { # goalname <- goals[1]
  cat(goalname, '\n')
  tmp_df <- scores_rgn %>%
    filter(goal == goalname) %>%
    filter(!is.na(status) & !is.na(obs_future_status))
  
  lfs_lm_goal <- lm(obs_future_status ~ status + trend + pressures + resilience, data = tmp_df) %>%
    broom::tidy() %>%
    mutate(goal  = goalname,
           n_obs = nrow(tmp_df))
  
  stat_tr_lm <- lm(obs_future_status ~ status + trend, data = tmp_df) %>%
    broom::tidy() %>%
    mutate(goal  = goalname,
           n_obs = nrow(tmp_df))
  
  lfs_lm_all <- bind_rows(lfs_lm_all, lfs_lm_goal)
  stat_tr_all <- bind_rows(stat_tr_all, stat_tr_lm)

}

## AO 
## FIS 
## CPP 
## CSS 
## CW 
## HAB 
## SPP 
## ICO 
## LSP 
## TR

lfs_lm_all <- lfs_lm_all %>%
  mutate(estimate = round(estimate, 4),
         std.error = round(std.error, 4),
         statistic = round(statistic, 4))
  
stat_tr_all <- stat_tr_all %>%
  mutate(estimate = round(estimate, 4),
         std.error = round(std.error, 4),
         statistic = round(statistic, 4))

write_csv(lfs_lm_all, 'int/lfs_lm_all.csv')
write_csv(stat_tr_all, 'int/stat_tr_all.csv')

DT::datatable(lfs_lm_all, caption = 'Likely Future Status model')

DT::datatable(stat_tr_all, caption = 'Status and trend only model')

0.4 Data layer year spans

Clipped to 1990 and later; some data layers go back farther but these will not typically inform scores except as reference points.

layer_targets <- read_csv(file.path(dir_calc, 'explore/int/layers_targets_master.csv')) %>%
  select(-target_element, -dimension) %>%
  distinct()

data_years <- read_csv(file.path(dir_calc, 'master/all_data_years.csv'))

# no_year_spans <- layer_targets %>% 
#   filter(!layer %in% data_years$layer_name) %>%
#   group_by(layer) %>%
#   summarize(targets = paste(target, collapse = ', '))
# 
# knitr::kable(no_year_spans) %>% paste(collapse = '')

year_spans <- data_years %>%
  full_join(layer_targets, by = c('layer_name' = 'layer')) %>%
  filter(target != 'spatial') %>%
  group_by(layer_name) %>%
  filter(year >= 1990) %>%
  summarize(year_min = min(year),
            year_max = max(year),
            targets = paste(unique(target) %>% sort(), collapse = ', ')) %>%
  ungroup() %>%
  arrange(layer_name) %>%
  mutate(layer_name = factor(layer_name, levels = rev(.$layer_name %>% unique), ordered = TRUE))

span_plot <- ggplot(year_spans, aes(x = layer_name, color = targets)) +
  ggtheme_plot(base_size = 8) +
  geom_linerange(aes(ymin = year_min, ymax = year_max), show.legend = FALSE) +
  labs(x = 'Layer name',
       y = 'Data year') +
  scale_color_manual(values = rep(brewer.pal(n = 8, name = 'Dark2'), 4)) +
  geom_text(aes(y = year_max, label = targets), color = 'grey20', size = 1.6, 
            vjust = 0, nudge_x = 0.1, hjust = 1, show.legend = FALSE) +
  coord_flip()

ggsave(file.path(dir_calc, 'explore/layers_data_years.png'), height = 6, width = 5)